We study a network of 500 coupled modified van der Pol oscillators. The value of a parameter associated with each oscillator is drawn from a normal distribution, giving a heterogeneous network. For strong enough coupling the oscillators all have the same period, and we consider periodic forcing of the network when it is in this state. By exploiting the correlations that quickly develop between the state of an oscillator and the value of its parameter we obtain an approximate low-dimensional description of the system in terms of the first few coefficients in a polynomial chaos expansion. Standard bifurcation analysis can then be performed on this low-dimensional system, and the results obtained from this predict very well the behaviour of the high-dimensional system for any set of realisations of the random parameter. Situations in which the method begins to fails are also discussed.